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ABSTRACT 



. Radiative cooling is an important ingredient in hydrodynamical models involving 

evolution of high temperature plasmas. Unfortunately, calculating an accurate cooling 
coefficient generally requires the solution of over a hundred differential equations to 
follow the ionization. We discuss here a simple 2-parameter approximation for the 
' cooling coefficient due to elements heavier than H and He, for the temperature range 

T = 10 4 — 10 8 K. Tests of the method show that it successfully tracks the ionization level 
in severe dynamical environments, and accurately approximates the non-equilibrium 
cooling coefficient of the trace elements, usually to within 10% in all cases for which 
cooling is actually important. The error is large only when the temperature is dropping 
so rapidly due to expansion that radiative cooling is negligible, but even in this situation, 
the ionization level is followed sufficiently accurately. The current approximation is fully 



o 

implemented in publicly available FORTRAN code. A second paper will discuss general 
approaches to approximation methods of this type, other realizations which could be 
even more accurate, and the potential for extension to calculations of non-equilibrium 
spectra. 

H- 

Subject headings: atomic data — radiation mechanisms:thermal — plasmas — hydro 
dynamics — supernova remnants — ISM:general 



1. Introduction 

In this paper we present a reliable but compact approximation to the high temperature cool- 
ing coefficient of heavy elements in diffuse plasmas under non-equilibrium conditions. 2 Such an 



1 Current address: Department of Physics, Stanford University, Stanford, CA 94305 

2 A set of FORTRAN routines to implement this approximation are provided at 
ftp://wisp5.physics.wisc.edu/pub/benjamin/BBC. Documentation is given in the file README. bbc. 
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approximation is essential for use with multi-dimensional hydrodynamics codes where the addi- 
tional burden of following the details of ionization evolution severely restricts the available spatial 
resolution. 

With only a few exceptions, large hydrodynamic models that incorporate radiative cooling 
characterize the cooling coefficient with a single parameter, the temperature, L(T e ), where the 
total emissivity per unit volume is A(T e ) = n e riHL(T e ) (ergs s" 1 cm~ 3 ) and n e and nn are the 
electron and hydrogen densities. These cooling functions are determined by assuming either that 
the ionization state at a given temperature is characterized by collisional equilibrium, or that all gas 
follows a particular pre-calculated ionization history (Shapiro & Moore 1976; Edgar Sz Chevalier 
1983; Sutherland & Dopita 1993). 

Because the cooling of a plasma depends on the ionization history of the constituent ions, 
there can actually be a large range in the value of L at a given T e , depending upon the details 
of the ionization evolution. We demonstrate, however, that cooling due to trace elements (those 
heavier than helium) can be approximated by following the evolution of just a single additional 
parameter, the mean charge on the trace ions, z = (Y,z, z Az z cz, z ) l^zAz- Here, Z is the element 
number, z is the ionic charge, Az is the (linear) abundance of element Z relative to hydrogen, 
and cz, z = nz, z /nz is the concentration of a given ion (the fraction of element Z with charge z). 
Because the abundance is dominated by oxygen, z ranges from zero to about nine. 

Our method generalizes from a cooling curve, L(T e ), to a cooling plane L(T e , z). It consequently 
requires another function that allows us to update the ionization level, D(T e ,z) = (1 / n e )(dz / dt) . 
With the mean charge as our ionization level indicator, D(T e , z) is the difference between the mean 
ionization and recombination functions. Thus D(T e ,z) = I(T e ,z) — R(T e ,z), where I(T e ,z) = 
(T, Z , Z A Z iz,z(T e ) c z ,z) l^zA z and R(T e ,z) = (T, Z , Z A Z r z , z (T e ) c z , z ) /^zA z , and where iz, z (T e ) 
and rz, z (T e ) are the ionization and recombination rate coefficients, respectively, at T e for stage z of 
element Z. The corresponding cooling coefficient is L(T e , z) = Y>z, z Az jz, z (T e ) c ZjZ , where jz, z {T e ) 
is the cooling coefficient per ion Z, z. 

To implement this approximation we search for a reasonable description of the nonequilibrium 
ionization concentrations, c ZtZ (T e , z), representative of those found at T e and z. This assumes 
that under actual nonequilibrium conditions, the distribution over ion states will depend less on 
the details of the past history than on how far the present mean charge differs from the collisional 
equilibrium value at the current temperature. In this paper, we consider the ionization distributions 
that arise in isothermal relaxation to equilibrium for each temperature T e , starting with nearly fully 
ionized or fully neutral gas. 

The atomic data used to develop the manifolds comes from Raymond & Smith (1977), with 
updates described in Raymond & Cox (1985), corrections in oscillator strengths in the cooling 
transitions of Li- like ions (Shull & Slavin 1994), and revised dielectronic recombination rates of 
Romanik (1988). The abundances are taken from Anders & Grevesse (1989). The effects of charge 
exchange are not included; errors introduced by neglecting this are smaller than the uncertainties 
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in the atomic data and elemental abundances used to generate the cooling curve. 

2. A Manifold of Isothermal Evolutions 

We have performed the manifold of non-equilibrium isothermal evolutions described below to 
form the basis for our approximation, producing tables of the cooling coefficient of trace elements 
L(T e ,z) and our ionization evolution functions I(T e ,z) and R(T e ,z). 

Gas is initialized with equilibrium ionization appropriate to some initial temperature T Q . Its 
temperature is then suddenly changed to T e and held fixed as the ionization evolves to equilib- 
rium. By choosing T Q to be low (e.g., 10 40 K), and performing the calculation for a dense set of 
temperatures T e between 10 4 and 10 80 K, the runs sample the full range of conditions possible 
in under-ionized gases. By repeating the whole set once again with very high T a (e.g. 10 80 K), 
conditions representing over-ionized gases are explored. Together these cases sample the full range 
of z. 

A particular moment in an evolution can be characterized by T Q , T e , and the fluence (/ = 
J n e dt). The ionization state is known, allowing straightforward evaluation of the mean charge z 
and then using z as the index of the time evolution. Evolution at constant temperature makes it 
easy to acquire the tables of L(T e , z), I(T e , z), R(T e , z), and D(T e , z) on a fixed grid of T e and z. 

The cooling function shown in Figure 1 shows a huge peak at T e ~ 10 6 K and z ~ 2 due 
to collisional excitation of low stages of ionization. It falls rapidly to the left due to Boltzmann 
factors, and gradually to the right due to decreasing excitation cross-sections. It falls at higher z as 
fewer bound electrons with low enough energy are available for collisional excitation. Above the 
equilibrium line the gas is over-ionized and collisional excitation is difficult, particularly at lower 
temperatures. At a given z, the mean ionization function, I, behaves much like the collisionally 
excited cooling, with vagaries at the transition from z = 1 to 2, and at z = 6 where the mean 
ionization rate is sensitive to the relative proportions of helium- and lithium-like oxygen. 

The structure of the mean recombination rate, R, is more gradual. In the upper left corner 
it is dominated by radiative recombination and the gradual increase with z and decrease with 
temperature are as expected. Deviations from that smooth pattern in the rest of the diagram are 
due to dielectronic recombination. In both I and R, there are distortions in the patterns just above 
the equilibrium line. These arise because, in relaxing to equilibrium, a recombining plasma goes 
through a considerable compaction of its ionization distribution. 

The principal feature of the rate function D(z, T e ) is that ionization toward equilibrium from 
below is much faster than recombination toward it from above. With increasing distance from 
equilibrium, the ionization rate increases dramatically, the recombination rate only gradually. 
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Fig. 1. — Contours of the trace element cooling coefficient, L(T e ,z) (erg cm 3 s _1 ) [top left panel], 
the mean ionization rate coefficient, I(T e ,z) (cm 3 s~ 4 ) [top right panel], the mean recombination 
rate coefficient, R(T e ,z) (cm 3 s _1 ) [bottom right panel], and the absolute value of the ionization 
evolution function, \D(z,T e )\ (cm 3 s -1 ) [bottom left panel], as a function of electron temperature 
and mean charge z. The dotted line in each panel shows the mean charge in collisional equilibrium as 
a function of gas temperature, z eq (T e ). The function \D(z, T e )\ (cm 3 s _1 ) passes through zero along 
this line. Above the equlibrium curve, the gas is recombining and D(z,T e ) is negative; below, it is 
ionizing and D(z,T e ) is positive. These functions were calculated for a grid of isothermal evolutions. 
Gas was either started with ionization fractions corresponding to T Q = I 4 K or T Q = 1 8 K. It then 
ionized up or recombined down to z eq (T e ) at fixed temperature T e . The resultant rate coefficients 
were calculated versus z for each T e for the range T e = 10 4 - I0 8 K. 

3. Testing With Cooling Blast Wave Evolutions 

To test the accuracy of our approximation for a wide range of situations, a representative set of 
cases were examined in which single parcels of gas in collisional equilibrium at 10 4 K were suddenly 
shock heated and then subjected to varying degrees of expansion. The test situation was modeled 
on the Sedov blast wave solution for an explosion of energy 10 51 ergs into a homogeneous medium 
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of particle density n . By varying the preshock density and the assumed distance of the parcel from 
the explosion site, we adjusted the post shock temperature and the timescales for depressurization 
versus radiative cooling. These scenarios test the two behaviors in which drastic departures from 
equilibrium occur: the rapid ionization of gas passing through a shock front and the rapid adiabatic 
cooling of an over-ionized plasma. 

Each parcel's evolution can be fully characterized by n Q and Ti- We have approximated 
p(t) = P2(t2/t) a , with t2 being the time the parcel is shocked, and a = 1.9 determining the rate of 
decompression. The temperature evolution is then 



where nu is the particle density of hydrogen (ionized and neutral) . 

We have carried out each simulation twice, first, solving the full ionization balance evolution 
exactly, then using our cooling approximation. In the latter, we solve the ionization evolution for 
hydrogen, helium, and z (four rate equations). The cooling is the sum of that from hydrogen, 
helium, and our tabulated cooling function. 

In order to test the widest possible range of conditions, we chose three values of T2, with 
three different values of n Q for each. The values of density were spaced to alter the initial ratio 
of depressurization to radiative cooling from approximately 0.4 at the high density end, to 6 and 
then 4000 at the lower densities. For the densest runs, the structure is very similar to that of a 
steady state radiative shock. For the next lower density cases, depressurization has strong effects 
but cooling is still sufficiently rapid that the density eventually rises rather than falling with time. 
At the lowest densities, the density and temperature fall adiabatically and the ionization level is 
soon frozen in. 

Figure 2 shows the exact trace element cooling coefficient evolutions versus temperature, the 
ratio of our approximate results to the exact ones, and the evolutions of mean charge, both exact 
and approximate. In each case, the cooling function begins very high just after the shock, when 
the gas is still briefly nearly neutral. At the two higher temperatures in the middle and right hand 
panels, the gas rapidly ionizes up and the cooling function drops precipitously to very close to 
the equilibrium value. This "ion flash" is so rapid that the integrated cooling during this period is 
negligible, as evidenced by the nearly vertical tracks. Thereafter, these tracks follow the equilibrium 
curve down to just below 10 6 K, after which their recombination cannot keep up with the cooling 
and the cooling rate falls below the equilibrium value. In the lowest temperature cases, significant 
cooling takes place before the cooling function drops through equilibrium, after which it plunges 
toward the nonequilibrium curve of the initially higher temperature cases. The differences between 
the radiation dominated and mixed evolutions are almost negligible. In the expansion dominated 
cases, however, things are quite different. The temperature drops so rapidly that the gas is soon 
over-ionized; recombination is slower than the temperature drop, and the ionization level declines 
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Fig. 2. — Comparison of cooling and ionization evolution for nine cases using the full ionization 
evolution and the approximation presented in this paper. The top row shows the evolution of 
the exact trace element cooling coefficient versus temperature. The left panel shows expansion- 
dominated evolution with (n , T 2 ) = (1 cur 3 , 3 x 10 7 K) (dashed), (3 x 1CT 4 cm" 3 , 3 x 10 6 K) (dash- 
dot), and (10~ 7 cm~ 3 , 3 x 10 5 K) (solid); the central panel shows "mixed" evolution with (n , T 2 ) = 
(2 x 10 4 cm' 3 , 3 x 10 7 K) (dashed), (10 cm" 3 , 3 x 10 6 K) (dash-dot), and (2 x 1(T 3 cm" 3 , 3 x 10 5 K) 
(solid); the right panel shows radiative cooling dominated evolution with (n ,T 2 ) = (10 6 cm~ 3 ,3 x 
10 7 K) (dashed), (300 cm" 3 , 3 x 10 6 K) (dash-dot), and (0.1 cm" 3 , 3 x 10 5 K) (solid). The time 
evolution for each curve is from right to left. The cooling curve for collisional ionization equilibrium 
is also shown (dotted). The middle set of panels shows the ratio of the approximation for trace 
element cooling to the exact value for each of the above cases. For "mixed" and "radiative" 
evolution, the highest temperature cases are sufficiently similar that the curves overlie each other. 
The lowest set of panels show the ionization evolution for these cases in the T-z plane. The solid 
curves show the exact evolution for the nine cases above, while the dashed lines show the slight 
departures that result by using the approximation we present. The dotted line in each panel shows 
the mean charge in collisional equilibrium as a function of gas temperature, z eq (T e ), as in Figure 1. 
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only gradually or is frozen in. 

The middle row of panels shows the ratio of approximate to exact cooling coefficients. In the 
two right hand panels, the approximation is good to within about 10% except below 10 42 K where 
it is consistently high by about 40%, and in a peak at 10 5,5 K where it can be as much as 20 to 
30% high. At this temperature, the approximate z evolutions of the two higher temperature cases 
in the two right lower panels of Figure 2 diverge by +0.4 charge states from the exact results. This 
is due in part to the excess cooling, but the recombination rate is almost certainly slightly too low 
as well. 

For the rapidly depressurizing cases, radiative cooling is inconsequential. Nevertheless, the 
cooling coefficient which differs substantially from the other cases, is still reasonably well fit, except 
for the intermediate temperature case for which z gets frozen at a value of 5.5. In that case, the 
value of z is well represented by the approximation, but the approximate cooling coefficient is too 
high. The freezing in of the ionization is extremely well matched, except very late in the highest 
initial temperature case. 

We have also compared the time evolution of density, temperature, mean charge, and the 
cooling coefficients for the exact calculation, our approximation, and the case where the trace 
elements are assumed to be in collisional ionization equilibrium. While the time evolution of 
temperature and density of the exact calculations and our approximation are in reasonably good 
agreement, the case with collisional ionization equilibrium shows significant differences in the time 
evolution of temperature. In addition, the ionization equilibrium assumption provides a poor 
estimate of the true ionization level of the gas, especially for the expansion dominated cases. 

4. Summary and Conclusions 

We provide a means by which large hydrocodes can include an accurate approximation to the 
radiative cooling coefficient, one far more responsive to the vagaries of dynamical environments 
than any single function of temperature. The method also follows the mean ionization level of the 
gas. Both the trace element cooling coefficient L(T e ,z), and the rate of change of z depend only 
on the identification of a manifold of representative ionization concentrations found at a particular 
temperature and mean charge. The ionization concentrations we used came from the isothermal 
relaxations of both highly over-ionized and highly under-ionized gas. In the test cases presented, 
both the cooling coefficient and mean charge evolution are quite accurately approximated. For all 
cases in which radiative cooling is significant, the approximation never errs by more than 30%. The 
error is usually less than 10 to 15%, well within the accuracy of the true cooling coefficient, which 
is limited by uncertainties in atomic data and abundances. The error in the cooling coefficient can 
be somewhat larger in examples with extreme amounts of depressurization, but never by enough 
to make the negligible radiative cooling appear to be significant. 

Our approximate cooling coefficients and charge evolution rates, however, showed patterns 
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in their modest errors, patterns which we believe we understand and can eliminate with future 
work. Should a potential user wait for these improvements before beginning to implement this 
method? The current model provides an excellent approximation for the cooling coefficient. If one is 
interested only in dynamics, it is certainly sufficient. Implementation of the next generation requires 
only swapping one set of tables for another. Our future work will also examine the possibility that 
the ion distributions can be used to provide absorption and emission spectra in post processing. It 
is possible to do this at the two parameter level, but higher order corrections will also be examined. 
The latter are expected to be two to four times more complex, and will be pursued only if spectral 
accuracy requires it. 

An additional consideration for gas with T ~ 10 4 K and z<2 is the effect of photoionization on 
the ionization balance and therefore the cooling. Our scheme can be easily modified to incorporate 
photoionization of hydrogen and helium. However, the effects of photoionization of trace elements 
would require a table of photoionization correction factors for L(T e , z) and D(T e , z). This too will 
be considered in future refinements. 

We would like to thank John Raymond for compiling and providing much of the atomic data 
that went into this work, and for a very careful reading of the manuscript. We would also like 
to thank NASA Astrophysical Theory grant NAG5-8417 for financial support of this work. And 
finally, we would like to acknowledge the valuable contributions of Leo Krzewina, Angela Klohs, 
Andrew Pawl, and Tim Freyer who all expended some effort on cracking this problem. 
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